{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8f40cb8f-9104-48ce-8ee6-fffa4671045c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "import pandas as pd\n",
    "from scipy.stats import spearmanr\n",
    "from statsmodels.sandbox.stats.multicomp import multipletests\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set_style('whitegrid')\n",
    "\n",
    "def p_adjust(pvalues, method='fdr_bh'):\n",
    "    res = multipletests(pvalues, method=method)\n",
    "    return np.array(res[1], dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3af6be6-d931-4378-9330-a064cbc0fb22",
   "metadata": {},
   "source": [
    "# Correlate nasal 16S with vaccine response\n",
    "\n",
    "##### Michael Shaffer\n",
    "##### 7/21/22\n",
    "##### Merck ESC, Sys bio group\n",
    "\n",
    "To look for associations between the nasal microbiome and vaccine response we have calculated correlations between the abundances of individual OTUs and the continuous titer measurements from 1 year of life."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1d18d0e-aa22-4f85-be4e-ed4c6f4a2516",
   "metadata": {},
   "source": [
    "## Read in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "663ff2b9-d055-41a8-80c7-70179b7c7872",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SubmissionType</th>\n",
       "      <th>SampleNumber</th>\n",
       "      <th>SampleIDValidation</th>\n",
       "      <th>DiversigenCheckInSampleName</th>\n",
       "      <th>ReplacesLowVolumeSampleID</th>\n",
       "      <th>BoxLocation</th>\n",
       "      <th>SampleType</th>\n",
       "      <th>SampleSource</th>\n",
       "      <th>SequencingType</th>\n",
       "      <th>BabyN</th>\n",
       "      <th>...</th>\n",
       "      <th>PCV ST9V_mmNorm</th>\n",
       "      <th>PCV ST14_mmNorm</th>\n",
       "      <th>PCV ST18C_mmNorm</th>\n",
       "      <th>PCV ST19A_mmNorm</th>\n",
       "      <th>PCV ST19F_mmNorm</th>\n",
       "      <th>PCV ST23F_mmNorm</th>\n",
       "      <th>median_mmNorm</th>\n",
       "      <th>median_mmNorm_DTAPHib</th>\n",
       "      <th>median_mmNorm_PCV</th>\n",
       "      <th>VR_group</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SampleID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>103_V5_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, A1</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>103</td>\n",
       "      <td>...</td>\n",
       "      <td>0.277753</td>\n",
       "      <td>0.146465</td>\n",
       "      <td>0.086356</td>\n",
       "      <td>0.035467</td>\n",
       "      <td>0.057453</td>\n",
       "      <td>0.300910</td>\n",
       "      <td>0.190854</td>\n",
       "      <td>0.326317</td>\n",
       "      <td>0.146465</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>108_V4_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, A9</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>108</td>\n",
       "      <td>...</td>\n",
       "      <td>0.871843</td>\n",
       "      <td>0.366194</td>\n",
       "      <td>0.164958</td>\n",
       "      <td>0.041692</td>\n",
       "      <td>0.063590</td>\n",
       "      <td>0.371471</td>\n",
       "      <td>0.168430</td>\n",
       "      <td>0.172291</td>\n",
       "      <td>0.164958</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>108_V5_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>7</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, B1</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>108</td>\n",
       "      <td>...</td>\n",
       "      <td>0.871843</td>\n",
       "      <td>0.366194</td>\n",
       "      <td>0.164958</td>\n",
       "      <td>0.041692</td>\n",
       "      <td>0.063590</td>\n",
       "      <td>0.371471</td>\n",
       "      <td>0.168430</td>\n",
       "      <td>0.172291</td>\n",
       "      <td>0.164958</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>117_V3_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>12</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, E1</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>117</td>\n",
       "      <td>...</td>\n",
       "      <td>0.082399</td>\n",
       "      <td>0.196600</td>\n",
       "      <td>0.114540</td>\n",
       "      <td>0.007407</td>\n",
       "      <td>0.003003</td>\n",
       "      <td>0.015946</td>\n",
       "      <td>0.114540</td>\n",
       "      <td>0.283004</td>\n",
       "      <td>0.066531</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>202_V5_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>14</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 1, E5</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>202</td>\n",
       "      <td>...</td>\n",
       "      <td>0.353074</td>\n",
       "      <td>0.061921</td>\n",
       "      <td>0.156257</td>\n",
       "      <td>0.175745</td>\n",
       "      <td>0.058983</td>\n",
       "      <td>0.154029</td>\n",
       "      <td>0.086954</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.154029</td>\n",
       "      <td>LVR</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 70 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "               SubmissionType  SampleNumber  SampleIDValidation  \\\n",
       "SampleID                                                          \n",
       "103_V5_NS_A1  Primary in Tube             1                 NaN   \n",
       "108_V4_NS_A1  Primary in Tube             6                 NaN   \n",
       "108_V5_NS_A1  Primary in Tube             7                 NaN   \n",
       "117_V3_NS_A1  Primary in Tube            12                 NaN   \n",
       "202_V5_NS_A1  Primary in Tube            14                 NaN   \n",
       "\n",
       "             DiversigenCheckInSampleName ReplacesLowVolumeSampleID  \\\n",
       "SampleID                                                             \n",
       "103_V5_NS_A1                         NaN                       NaN   \n",
       "108_V4_NS_A1                         NaN                       NaN   \n",
       "108_V5_NS_A1                         NaN                       NaN   \n",
       "117_V3_NS_A1                         NaN                       NaN   \n",
       "202_V5_NS_A1                         NaN                       NaN   \n",
       "\n",
       "             BoxLocation  SampleType  SampleSource SequencingType  BabyN  ...  \\\n",
       "SampleID                                                                  ...   \n",
       "103_V5_NS_A1   Box 1, A1  Nasal Swab  Human Infant            16S    103  ...   \n",
       "108_V4_NS_A1   Box 1, A9  Nasal Swab  Human Infant            16S    108  ...   \n",
       "108_V5_NS_A1   Box 1, B1  Nasal Swab  Human Infant            16S    108  ...   \n",
       "117_V3_NS_A1   Box 1, E1  Nasal Swab  Human Infant            16S    117  ...   \n",
       "202_V5_NS_A1   Box 1, E5  Nasal Swab  Human Infant            16S    202  ...   \n",
       "\n",
       "              PCV ST9V_mmNorm PCV ST14_mmNorm  PCV ST18C_mmNorm  \\\n",
       "SampleID                                                          \n",
       "103_V5_NS_A1         0.277753        0.146465          0.086356   \n",
       "108_V4_NS_A1         0.871843        0.366194          0.164958   \n",
       "108_V5_NS_A1         0.871843        0.366194          0.164958   \n",
       "117_V3_NS_A1         0.082399        0.196600          0.114540   \n",
       "202_V5_NS_A1         0.353074        0.061921          0.156257   \n",
       "\n",
       "             PCV ST19A_mmNorm  PCV ST19F_mmNorm PCV ST23F_mmNorm  \\\n",
       "SampleID                                                           \n",
       "103_V5_NS_A1         0.035467          0.057453         0.300910   \n",
       "108_V4_NS_A1         0.041692          0.063590         0.371471   \n",
       "108_V5_NS_A1         0.041692          0.063590         0.371471   \n",
       "117_V3_NS_A1         0.007407          0.003003         0.015946   \n",
       "202_V5_NS_A1         0.175745          0.058983         0.154029   \n",
       "\n",
       "              median_mmNorm median_mmNorm_DTAPHib median_mmNorm_PCV  VR_group  \n",
       "SampleID                                                                       \n",
       "103_V5_NS_A1       0.190854              0.326317          0.146465       NVR  \n",
       "108_V4_NS_A1       0.168430              0.172291          0.164958       NVR  \n",
       "108_V5_NS_A1       0.168430              0.172291          0.164958       NVR  \n",
       "117_V3_NS_A1       0.114540              0.283004          0.066531       NVR  \n",
       "202_V5_NS_A1       0.086954              0.000000          0.154029       LVR  \n",
       "\n",
       "[5 rows x 70 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meta = pd.read_csv('../../data/metadata/nasal/nasal_metadata.csv', index_col='SampleID')\n",
    "meta['age_at_collection'] = (pd.to_datetime(meta['CollectionDate']) - pd.to_datetime(meta['DOB'])).dt.days\n",
    "meta = pd.concat([meta,\n",
    "                  pd.read_csv('../../data/metadata/nasal/nasal_abx_usage.csv', index_col='SampleID'),\n",
    "                  pd.read_csv('../../data/metadata/nasal/nasal_titers_yr2.csv', index_col='SampleID')],\n",
    "                 axis=1)\n",
    "meta = meta.loc[~pd.isna(meta['median_mmNorm'])]\n",
    "meta.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "07db70cc-f930-477f-ab3f-27777f8af4f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>101_S1_NS_A1</th>\n",
       "      <th>101_V3_NS_A1</th>\n",
       "      <th>101_V5_NS_A1</th>\n",
       "      <th>102_V1_NS_A1</th>\n",
       "      <th>102_V3_NS_A1</th>\n",
       "      <th>102_V5_NS_A1</th>\n",
       "      <th>102_V6_NS_A1</th>\n",
       "      <th>103_S1_NS_A1</th>\n",
       "      <th>103_S3_NS_A1</th>\n",
       "      <th>103_V10_NS_A1</th>\n",
       "      <th>...</th>\n",
       "      <th>MSA2002_5A</th>\n",
       "      <th>MSA2002_5B</th>\n",
       "      <th>MSA2002_6A</th>\n",
       "      <th>MSA2002_6B</th>\n",
       "      <th>MSA2002_7A</th>\n",
       "      <th>MSA2002_7B</th>\n",
       "      <th>MSA2002_8A</th>\n",
       "      <th>MSA2002_8B</th>\n",
       "      <th>MSA2002_9A</th>\n",
       "      <th>MSA2002_9B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0001</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1593</td>\n",
       "      <td>7320</td>\n",
       "      <td>606</td>\n",
       "      <td>...</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0002</th>\n",
       "      <td>5845</td>\n",
       "      <td>9876</td>\n",
       "      <td>692</td>\n",
       "      <td>557</td>\n",
       "      <td>783</td>\n",
       "      <td>509</td>\n",
       "      <td>6047</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>114</td>\n",
       "      <td>126</td>\n",
       "      <td>104</td>\n",
       "      <td>115</td>\n",
       "      <td>119</td>\n",
       "      <td>111</td>\n",
       "      <td>168</td>\n",
       "      <td>147</td>\n",
       "      <td>83</td>\n",
       "      <td>103</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0003</th>\n",
       "      <td>117</td>\n",
       "      <td>0</td>\n",
       "      <td>879</td>\n",
       "      <td>4392</td>\n",
       "      <td>1428</td>\n",
       "      <td>528</td>\n",
       "      <td>87</td>\n",
       "      <td>877</td>\n",
       "      <td>2642</td>\n",
       "      <td>1498</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0004</th>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>1104</td>\n",
       "      <td>1</td>\n",
       "      <td>6133</td>\n",
       "      <td>475</td>\n",
       "      <td>1</td>\n",
       "      <td>109</td>\n",
       "      <td>2</td>\n",
       "      <td>14</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0005</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4173</td>\n",
       "      <td>24</td>\n",
       "      <td>3121</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 980 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         101_S1_NS_A1  101_V3_NS_A1  101_V5_NS_A1  102_V1_NS_A1  102_V3_NS_A1  \\\n",
       "Otu0001             1             0             0             1             2   \n",
       "Otu0002          5845          9876           692           557           783   \n",
       "Otu0003           117             0           879          4392          1428   \n",
       "Otu0004             9             1          1104             1          6133   \n",
       "Otu0005             0             0             0             0             0   \n",
       "\n",
       "         102_V5_NS_A1  102_V6_NS_A1  103_S1_NS_A1  103_S3_NS_A1  \\\n",
       "Otu0001             0             2          1593          7320   \n",
       "Otu0002           509          6047             1             0   \n",
       "Otu0003           528            87           877          2642   \n",
       "Otu0004           475             1           109             2   \n",
       "Otu0005             1             0          4173            24   \n",
       "\n",
       "         103_V10_NS_A1  ...  MSA2002_5A  MSA2002_5B  MSA2002_6A  MSA2002_6B  \\\n",
       "Otu0001            606  ...           3           4           1           1   \n",
       "Otu0002              3  ...         114         126         104         115   \n",
       "Otu0003           1498  ...           0           0           0           0   \n",
       "Otu0004             14  ...           0           0           1           0   \n",
       "Otu0005           3121  ...           0           0           0           0   \n",
       "\n",
       "         MSA2002_7A  MSA2002_7B  MSA2002_8A  MSA2002_8B  MSA2002_9A  \\\n",
       "Otu0001           0           2           1           1           1   \n",
       "Otu0002         119         111         168         147          83   \n",
       "Otu0003           0           0           0           0           0   \n",
       "Otu0004           0           0           0           0           0   \n",
       "Otu0005           0           0           0           0           0   \n",
       "\n",
       "         MSA2002_9B  \n",
       "Otu0001           0  \n",
       "Otu0002         103  \n",
       "Otu0003           0  \n",
       "Otu0004           0  \n",
       "Otu0005           0  \n",
       "\n",
       "[5 rows x 980 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts = pd.read_csv('../../data/nasal/otu_table.gt10_rar10K.tsv', sep='\\t', index_col=0).transpose()\n",
    "counts.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "be1270b2-e08a-4060-a5e3-3008288626f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(645, 70)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/6t/1w2t3qmd1rx81mfw9sq_tfpr0000gn/T/ipykernel_12545/52385233.py:2: FutureWarning: Passing a set as an indexer is deprecated and will raise in a future version. Use a list instead.\n",
      "  meta = meta.loc[in_both].sort_values(['BabyN', 'age_at_collection'])\n"
     ]
    }
   ],
   "source": [
    "in_both = set(meta.index) & set(counts.columns)\n",
    "meta = meta.loc[in_both].sort_values(['BabyN', 'age_at_collection'])\n",
    "print(meta.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "19ccfe3f-9562-46f7-8adf-91892eadfe82",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta_v5 = meta.query(\"VisitCode == 'V5'\")\n",
    "counts_v5 = counts[meta_v5.index]\n",
    "counts_v5 = counts_v5.loc[(counts_v5 > 0).sum(axis=1) > counts_v5.shape[1]*.2]\n",
    "\n",
    "meta_v6 = meta.query(\"VisitCode == 'V6'\")\n",
    "counts_v6 = counts[meta_v6.index]\n",
    "counts_v6 = counts_v6.loc[(counts_v6 > 0).sum(axis=1) > counts_v6.shape[1]*.2]\n",
    "\n",
    "meta_v7 = meta.query(\"VisitCode == 'V7'\")\n",
    "counts_v7 = counts[meta_v7.index]\n",
    "counts_v7 = counts_v7.loc[(counts_v7 > 0).sum(axis=1) > counts_v7.shape[1]*.2]\n",
    "\n",
    "meta_v9 = meta.query(\"VisitCode == 'V9'\")\n",
    "counts_v9 = counts[meta_v9.index]\n",
    "counts_v9 = counts_v9.loc[(counts_v9 > 0).sum(axis=1) > counts_v9.shape[1]*.2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7c1b3fe-5a38-4624-b29d-d1fb437276c8",
   "metadata": {},
   "source": [
    "There are some samples which have DTAPHib titers but do not have PCV titers. So here we refilter everything to remove samples that do not have PCV titers for analysis of PCV titers specifically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5c6245d0-c5fa-4693-bab7-d570d30938ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta_PCV = meta.loc[~pd.isna(meta['median_mmNorm_PCV'])]\n",
    "meta_PCV.head()\n",
    "\n",
    "meta_PCV_v5 = meta_PCV.query(\"VisitCode == 'V5'\")\n",
    "counts_PCV_v5 = counts[meta_PCV_v5.index]\n",
    "counts_PCV_v5 = counts_PCV_v5.loc[(counts_PCV_v5 > 0).sum(axis=1) > counts_PCV_v5.shape[1]*.2]\n",
    "\n",
    "meta_PCV_v6 = meta_PCV.query(\"VisitCode == 'V6'\")\n",
    "counts_PCV_v6 = counts[meta_PCV_v6.index]\n",
    "counts_PCV_v6 = counts_PCV_v6.loc[(counts_PCV_v6 > 0).sum(axis=1) > counts_PCV_v6.shape[1]*.2]\n",
    "\n",
    "meta_PCV_v7 = meta_PCV.query(\"VisitCode == 'V7'\")\n",
    "counts_PCV_v7 = counts[meta_PCV_v7.index]\n",
    "counts_PCV_v7 = counts_PCV_v7.loc[(counts_PCV_v7 > 0).sum(axis=1) > counts_PCV_v7.shape[1]*.2]\n",
    "\n",
    "meta_PCV_v9 = meta_PCV.query(\"VisitCode == 'V9'\")\n",
    "counts_PCV_v9 = counts[meta_PCV_v9.index]\n",
    "counts_PCV_v9 = counts_PCV_v9.loc[(counts_PCV_v9 > 0).sum(axis=1) > counts_PCV_v9.shape[1]*.2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03b6ca75-7f7f-4ffb-a575-5d399ab746e5",
   "metadata": {},
   "source": [
    "## Correlations with median titer values\n",
    "\n",
    "We will use Spearman's R as our correlation metric and use OTU abundances from the 2 month (V5), 4 month (V6), 6 month (V7) and 1 year (V9) time points. 2, 4 and 6 months are when vaccinations are given and 1 year is when titers were measured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "07c9fcf2-97a8-47a9-833d-e3eee803110d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0037</th>\n",
       "      <td>-0.316588</td>\n",
       "      <td>0.026674</td>\n",
       "      <td>0.368315</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0006</th>\n",
       "      <td>0.309796</td>\n",
       "      <td>0.030297</td>\n",
       "      <td>0.368315</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0035</th>\n",
       "      <td>-0.307117</td>\n",
       "      <td>0.031833</td>\n",
       "      <td>0.368315</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0079</th>\n",
       "      <td>-0.300463</td>\n",
       "      <td>0.035933</td>\n",
       "      <td>0.368315</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0009</th>\n",
       "      <td>-0.280625</td>\n",
       "      <td>0.050807</td>\n",
       "      <td>0.405400</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0037 -0.316588  0.026674  0.368315\n",
       "Otu0006  0.309796  0.030297  0.368315\n",
       "Otu0035 -0.307117  0.031833  0.368315\n",
       "Otu0079 -0.300463  0.035933  0.368315\n",
       "Otu0009 -0.280625  0.050807  0.405400"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v5_correlations = counts_v5.transpose().apply(spearmanr, b=meta_v5['median_mmNorm']).transpose()\n",
    "v5_correlations.columns = ['rho', 'p_value']\n",
    "v5_correlations['p_adj'] = p_adjust(v5_correlations['p_value'])\n",
    "v5_correlations = v5_correlations.sort_values('p_value')\n",
    "v5_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5d4719b9-ee9c-4fae-9806-cf6641bd571c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0009</th>\n",
       "      <td>-0.284144</td>\n",
       "      <td>0.052911</td>\n",
       "      <td>0.831504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0002</th>\n",
       "      <td>-0.260329</td>\n",
       "      <td>0.077181</td>\n",
       "      <td>0.831504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0050</th>\n",
       "      <td>-0.254280</td>\n",
       "      <td>0.084565</td>\n",
       "      <td>0.831504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0015</th>\n",
       "      <td>-0.247202</td>\n",
       "      <td>0.093894</td>\n",
       "      <td>0.831504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0011</th>\n",
       "      <td>-0.241853</td>\n",
       "      <td>0.101457</td>\n",
       "      <td>0.831504</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0009 -0.284144  0.052911  0.831504\n",
       "Otu0002 -0.260329  0.077181  0.831504\n",
       "Otu0050 -0.254280  0.084565  0.831504\n",
       "Otu0015 -0.247202  0.093894  0.831504\n",
       "Otu0011 -0.241853  0.101457  0.831504"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v6_correlations = counts_v6.transpose().apply(spearmanr, b=meta_v6['median_mmNorm']).transpose()\n",
    "v6_correlations.columns = ['rho', 'p_value']\n",
    "v6_correlations['p_adj'] = p_adjust(v6_correlations['p_value'])\n",
    "v6_correlations = v6_correlations.sort_values('p_value')\n",
    "v6_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "dab2c6b9-30ad-4be3-b858-2114a8308d00",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0049</th>\n",
       "      <td>0.325078</td>\n",
       "      <td>0.021253</td>\n",
       "      <td>0.73654</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0068</th>\n",
       "      <td>-0.288976</td>\n",
       "      <td>0.041817</td>\n",
       "      <td>0.73654</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0006</th>\n",
       "      <td>-0.278456</td>\n",
       "      <td>0.050219</td>\n",
       "      <td>0.73654</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0039</th>\n",
       "      <td>-0.255738</td>\n",
       "      <td>0.073046</td>\n",
       "      <td>0.77184</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0004</th>\n",
       "      <td>-0.217258</td>\n",
       "      <td>0.129632</td>\n",
       "      <td>0.77184</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value    p_adj\n",
       "Otu0049  0.325078  0.021253  0.73654\n",
       "Otu0068 -0.288976  0.041817  0.73654\n",
       "Otu0006 -0.278456  0.050219  0.73654\n",
       "Otu0039 -0.255738  0.073046  0.77184\n",
       "Otu0004 -0.217258  0.129632  0.77184"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v7_correlations = counts_v7.transpose().apply(spearmanr, b=meta_v7['median_mmNorm']).transpose()\n",
    "v7_correlations.columns = ['rho', 'p_value']\n",
    "v7_correlations['p_adj'] = p_adjust(v7_correlations['p_value'])\n",
    "v7_correlations = v7_correlations.sort_values('p_value')\n",
    "v7_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "42e16d24-3221-4620-ad29-da3021ed7d0c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0067</th>\n",
       "      <td>-0.252197</td>\n",
       "      <td>0.065804</td>\n",
       "      <td>0.919694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0123</th>\n",
       "      <td>0.212457</td>\n",
       "      <td>0.122983</td>\n",
       "      <td>0.919694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0040</th>\n",
       "      <td>0.208318</td>\n",
       "      <td>0.130627</td>\n",
       "      <td>0.919694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0090</th>\n",
       "      <td>0.202293</td>\n",
       "      <td>0.142383</td>\n",
       "      <td>0.919694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0014</th>\n",
       "      <td>0.199677</td>\n",
       "      <td>0.147727</td>\n",
       "      <td>0.919694</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0067 -0.252197  0.065804  0.919694\n",
       "Otu0123  0.212457  0.122983  0.919694\n",
       "Otu0040  0.208318  0.130627  0.919694\n",
       "Otu0090  0.202293  0.142383  0.919694\n",
       "Otu0014  0.199677  0.147727  0.919694"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v9_correlations = counts_v9.transpose().apply(spearmanr, b=meta_v9['median_mmNorm']).transpose()\n",
    "v9_correlations.columns = ['rho', 'p_value']\n",
    "v9_correlations['p_adj'] = p_adjust(v9_correlations['p_value'])\n",
    "v9_correlations = v9_correlations.sort_values('p_value')\n",
    "v9_correlations.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "967bc78d-674e-47de-a92f-902130f50e84",
   "metadata": {},
   "source": [
    "Across all test nothing is significant after multiple test correction. OTU 140 is significant raw at V6 and V7 but not strong enough to dig into further. The most raw significant results are at 1 year. Another potential indicator that titer is effected by factors present when it is measured?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "895a3234-112b-4ec5-8d2e-52bbb6d57f31",
   "metadata": {},
   "source": [
    "## Correlations with median titer group values\n",
    "\n",
    "Now we will split the titers into DTAPHib and PCV and test for significant correlations separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dfaa593b-4625-4c9a-935e-501ea409e857",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0039</th>\n",
       "      <td>-0.279394</td>\n",
       "      <td>0.051874</td>\n",
       "      <td>0.987112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0029</th>\n",
       "      <td>0.257390</td>\n",
       "      <td>0.074191</td>\n",
       "      <td>0.987112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0068</th>\n",
       "      <td>-0.219554</td>\n",
       "      <td>0.129578</td>\n",
       "      <td>0.987112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0001</th>\n",
       "      <td>0.213858</td>\n",
       "      <td>0.140081</td>\n",
       "      <td>0.987112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0025</th>\n",
       "      <td>-0.198642</td>\n",
       "      <td>0.171226</td>\n",
       "      <td>0.987112</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0039 -0.279394  0.051874  0.987112\n",
       "Otu0029  0.257390  0.074191  0.987112\n",
       "Otu0068 -0.219554  0.129578  0.987112\n",
       "Otu0001  0.213858  0.140081  0.987112\n",
       "Otu0025 -0.198642  0.171226  0.987112"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v5_DTAPHib_correlations = counts_v5.transpose().apply(spearmanr, b=meta_v5['median_mmNorm_DTAPHib']).transpose()\n",
    "v5_DTAPHib_correlations.columns = ['rho', 'p_value']\n",
    "v5_DTAPHib_correlations['p_adj'] = p_adjust(v5_DTAPHib_correlations['p_value'])\n",
    "v5_DTAPHib_correlations = v5_DTAPHib_correlations.sort_values('p_value')\n",
    "v5_DTAPHib_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2a6ff897-fec6-4a82-a4fc-e2d4de36c2e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0009</th>\n",
       "      <td>-0.381482</td>\n",
       "      <td>0.006840</td>\n",
       "      <td>0.180739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0079</th>\n",
       "      <td>-0.370310</td>\n",
       "      <td>0.008817</td>\n",
       "      <td>0.180739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0035</th>\n",
       "      <td>-0.316148</td>\n",
       "      <td>0.026897</td>\n",
       "      <td>0.208397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0021</th>\n",
       "      <td>-0.308617</td>\n",
       "      <td>0.030965</td>\n",
       "      <td>0.208397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0008</th>\n",
       "      <td>-0.305447</td>\n",
       "      <td>0.032824</td>\n",
       "      <td>0.208397</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0009 -0.381482  0.006840  0.180739\n",
       "Otu0079 -0.370310  0.008817  0.180739\n",
       "Otu0035 -0.316148  0.026897  0.208397\n",
       "Otu0021 -0.308617  0.030965  0.208397\n",
       "Otu0008 -0.305447  0.032824  0.208397"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v5_PCV_correlations = counts_PCV_v5.transpose().apply(spearmanr, b=meta_PCV_v5['median_mmNorm_PCV']).transpose()\n",
    "v5_PCV_correlations.columns = ['rho', 'p_value']\n",
    "v5_PCV_correlations['p_adj'] = p_adjust(v5_PCV_correlations['p_value'])\n",
    "v5_PCV_correlations = v5_PCV_correlations.sort_values('p_value')\n",
    "v5_PCV_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9f5e043a-fafc-402a-a8bc-de132b2a0821",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0002</th>\n",
       "      <td>-0.342277</td>\n",
       "      <td>0.018524</td>\n",
       "      <td>0.77051</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0079</th>\n",
       "      <td>-0.272942</td>\n",
       "      <td>0.063421</td>\n",
       "      <td>0.77051</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0063</th>\n",
       "      <td>0.232247</td>\n",
       "      <td>0.116206</td>\n",
       "      <td>0.77051</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0001</th>\n",
       "      <td>0.224774</td>\n",
       "      <td>0.128766</td>\n",
       "      <td>0.77051</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0026</th>\n",
       "      <td>-0.220130</td>\n",
       "      <td>0.137066</td>\n",
       "      <td>0.77051</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value    p_adj\n",
       "Otu0002 -0.342277  0.018524  0.77051\n",
       "Otu0079 -0.272942  0.063421  0.77051\n",
       "Otu0063  0.232247  0.116206  0.77051\n",
       "Otu0001  0.224774  0.128766  0.77051\n",
       "Otu0026 -0.220130  0.137066  0.77051"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v6_DTAPHib_correlations = counts_v6.transpose().apply(spearmanr, b=meta_v6['median_mmNorm_DTAPHib']).transpose()\n",
    "v6_DTAPHib_correlations.columns = ['rho', 'p_value']\n",
    "v6_DTAPHib_correlations['p_adj'] = p_adjust(v6_DTAPHib_correlations['p_value'])\n",
    "v6_DTAPHib_correlations = v6_DTAPHib_correlations.sort_values('p_value')\n",
    "v6_DTAPHib_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "72ccf352-6156-4ce1-9831-76cd256e3d9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0009</th>\n",
       "      <td>-0.320176</td>\n",
       "      <td>0.028234</td>\n",
       "      <td>0.85747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0021</th>\n",
       "      <td>-0.282048</td>\n",
       "      <td>0.054763</td>\n",
       "      <td>0.85747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0008</th>\n",
       "      <td>-0.271260</td>\n",
       "      <td>0.065134</td>\n",
       "      <td>0.85747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0011</th>\n",
       "      <td>-0.252378</td>\n",
       "      <td>0.086997</td>\n",
       "      <td>0.85747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0035</th>\n",
       "      <td>-0.246389</td>\n",
       "      <td>0.095015</td>\n",
       "      <td>0.85747</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value    p_adj\n",
       "Otu0009 -0.320176  0.028234  0.85747\n",
       "Otu0021 -0.282048  0.054763  0.85747\n",
       "Otu0008 -0.271260  0.065134  0.85747\n",
       "Otu0011 -0.252378  0.086997  0.85747\n",
       "Otu0035 -0.246389  0.095015  0.85747"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v6_PCV_correlations = counts_PCV_v6.transpose().apply(spearmanr, b=meta_PCV_v6['median_mmNorm_PCV']).transpose()\n",
    "v6_PCV_correlations.columns = ['rho', 'p_value']\n",
    "v6_PCV_correlations['p_adj'] = p_adjust(v6_PCV_correlations['p_value'])\n",
    "v6_PCV_correlations = v6_PCV_correlations.sort_values('p_value')\n",
    "v6_PCV_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9af93a96-e47c-4fc9-86d1-abb5e2db5dd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0044</th>\n",
       "      <td>0.266439</td>\n",
       "      <td>0.061439</td>\n",
       "      <td>0.997358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0034</th>\n",
       "      <td>-0.253045</td>\n",
       "      <td>0.076225</td>\n",
       "      <td>0.997358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0072</th>\n",
       "      <td>-0.234675</td>\n",
       "      <td>0.100916</td>\n",
       "      <td>0.997358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0036</th>\n",
       "      <td>-0.220739</td>\n",
       "      <td>0.123451</td>\n",
       "      <td>0.997358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0050</th>\n",
       "      <td>-0.208518</td>\n",
       "      <td>0.146173</td>\n",
       "      <td>0.997358</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0044  0.266439  0.061439  0.997358\n",
       "Otu0034 -0.253045  0.076225  0.997358\n",
       "Otu0072 -0.234675  0.100916  0.997358\n",
       "Otu0036 -0.220739  0.123451  0.997358\n",
       "Otu0050 -0.208518  0.146173  0.997358"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v7_DTAPHib_correlations = counts_v7.transpose().apply(spearmanr, b=meta_v7['median_mmNorm_DTAPHib']).transpose()\n",
    "v7_DTAPHib_correlations.columns = ['rho', 'p_value']\n",
    "v7_DTAPHib_correlations['p_adj'] = p_adjust(v7_DTAPHib_correlations['p_value'])\n",
    "v7_DTAPHib_correlations = v7_DTAPHib_correlations.sort_values('p_value')\n",
    "v7_DTAPHib_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "db30229f-5604-476f-a92c-960aba714300",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0049</th>\n",
       "      <td>0.392223</td>\n",
       "      <td>0.004845</td>\n",
       "      <td>0.210327</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0006</th>\n",
       "      <td>-0.357307</td>\n",
       "      <td>0.010855</td>\n",
       "      <td>0.210327</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0004</th>\n",
       "      <td>-0.344334</td>\n",
       "      <td>0.014340</td>\n",
       "      <td>0.210327</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0002</th>\n",
       "      <td>-0.232003</td>\n",
       "      <td>0.104971</td>\n",
       "      <td>0.740988</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0039</th>\n",
       "      <td>-0.228355</td>\n",
       "      <td>0.110705</td>\n",
       "      <td>0.740988</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0049  0.392223  0.004845  0.210327\n",
       "Otu0006 -0.357307  0.010855  0.210327\n",
       "Otu0004 -0.344334  0.014340  0.210327\n",
       "Otu0002 -0.232003  0.104971  0.740988\n",
       "Otu0039 -0.228355  0.110705  0.740988"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "v7_PCV_correlations = counts_PCV_v7.transpose().apply(spearmanr, b=meta_PCV_v7['median_mmNorm_PCV']).transpose()\n",
    "v7_PCV_correlations.columns = ['rho', 'p_value']\n",
    "v7_PCV_correlations['p_adj'] = p_adjust(v7_PCV_correlations['p_value'])\n",
    "v7_PCV_correlations = v7_PCV_correlations.sort_values('p_value')\n",
    "v7_PCV_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1ec6b34b-73cf-4121-90ae-3bdd5b157f56",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0093</th>\n",
       "      <td>-0.346945</td>\n",
       "      <td>0.010163</td>\n",
       "      <td>0.741886</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0008</th>\n",
       "      <td>-0.287901</td>\n",
       "      <td>0.034768</td>\n",
       "      <td>0.968702</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0002</th>\n",
       "      <td>-0.274830</td>\n",
       "      <td>0.044301</td>\n",
       "      <td>0.968702</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0095</th>\n",
       "      <td>-0.240592</td>\n",
       "      <td>0.079696</td>\n",
       "      <td>0.968702</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0003</th>\n",
       "      <td>0.211718</td>\n",
       "      <td>0.124323</td>\n",
       "      <td>0.968702</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0093 -0.346945  0.010163  0.741886\n",
       "Otu0008 -0.287901  0.034768  0.968702\n",
       "Otu0002 -0.274830  0.044301  0.968702\n",
       "Otu0095 -0.240592  0.079696  0.968702\n",
       "Otu0003  0.211718  0.124323  0.968702"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v9_DTAPHib_correlations = counts_v9.transpose().apply(spearmanr, b=meta_v9['median_mmNorm_DTAPHib']).transpose()\n",
    "v9_DTAPHib_correlations.columns = ['rho', 'p_value']\n",
    "v9_DTAPHib_correlations['p_adj'] = p_adjust(v9_DTAPHib_correlations['p_value'])\n",
    "v9_DTAPHib_correlations = v9_DTAPHib_correlations.sort_values('p_value')\n",
    "v9_DTAPHib_correlations.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7f2cd382-5979-429c-af3a-c6eb08d26b56",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rho</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Otu0014</th>\n",
       "      <td>0.356955</td>\n",
       "      <td>0.008059</td>\n",
       "      <td>0.588275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0138</th>\n",
       "      <td>0.304657</td>\n",
       "      <td>0.025096</td>\n",
       "      <td>0.916009</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0090</th>\n",
       "      <td>0.281898</td>\n",
       "      <td>0.038910</td>\n",
       "      <td>0.946810</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0067</th>\n",
       "      <td>-0.264554</td>\n",
       "      <td>0.053214</td>\n",
       "      <td>0.954329</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Otu0006</th>\n",
       "      <td>-0.235970</td>\n",
       "      <td>0.085833</td>\n",
       "      <td>0.954329</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              rho   p_value     p_adj\n",
       "Otu0014  0.356955  0.008059  0.588275\n",
       "Otu0138  0.304657  0.025096  0.916009\n",
       "Otu0090  0.281898  0.038910  0.946810\n",
       "Otu0067 -0.264554  0.053214  0.954329\n",
       "Otu0006 -0.235970  0.085833  0.954329"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v9_PCV_correlations = counts_PCV_v9.transpose().apply(spearmanr, b=meta_PCV_v9['median_mmNorm_PCV']).transpose()\n",
    "v9_PCV_correlations.columns = ['rho', 'p_value']\n",
    "v9_PCV_correlations['p_adj'] = p_adjust(v9_PCV_correlations['p_value'])\n",
    "v9_PCV_correlations = v9_PCV_correlations.sort_values('p_value')\n",
    "v9_PCV_correlations.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f36b9cd-bad9-4983-941e-3b13a77e5df3",
   "metadata": {},
   "source": [
    "Still nothing significant after multiple testing correction. Still most raw significance at 1 year (V9) but nothing consistent enough to get excited about despite the low p-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "512c00ba-d5ff-421b-a9f8-2779994ebc29",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
